On the optimal calculation of the pair correlation function for an orthorombic system 
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We present a new computational method to calculate arbitrary pair correlation functions of an 
orthorombic system in the most efficient way. The algorithm is demonstrated by the calculation of 
the radial distribution function of shock compressed liquid hydrogen. 



Determining the equilibrium properties of a classical 
many-body system is rather challenging due to the large 
variety of behavior, which results from the interactions 
between the particles. The atomic structure of a system 
is best characterized by a set of correlation functions 
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where D(V) is the spatial domain defined by the physi- 
cal container, N the number of particles, V the volume, 
(3 = 1/ksT the inverse temperature and Un(Ri, ■■■,R n ) 
the potential energy function of the system, while n is the 
number of fixed particles. The simplest is the so called 
pair-correlation function (PCF) g^(Ri,R 2 ) = g(r), 
which is of particular importance, as it is related to the 
structure factor S(q) = 1 + p § D ry) dr e~ lqr {g(r) — 1] by 
the Fourier transform. For a spatially isotropic system, 
such as a liquid, the PCF does only depend on the ab- 
solute value of the relative distance |r| = r and can be 
simplified to the radial distribution function (RDF) g(r), 
which gives the probability of finding a pair of atoms a 
distance r apart, relative to the probability of a ideal 
gas at the same density. Like before, the Fourier trans- 
form of g(r) corresponds to the isotropic structure factor 
5(|q|) = 1 + p f D , v ^ dr e~ lqr [g(r) — 1], which can be ex- 
perimentally measured using neutron scattering or X-ray 
diffraction techniques [TJ [2] . In addition, it is not only 
possible to calculate the potential of mean force, but as- 
suming that Un{Ri, ■■■,R n ) is a pair-potential, also nu- 
merous important thermodynamic properties such as the 
potential energy and the pressure, just to name a few. 

In order to compute the RDF in a periodic cubic cell 
of length L, it is customary to restrict oneself to consider 
only those atoms that lie within the sphere of radius L/2 
and discard all the rest. In this most favorable case, the 
simulation cell needs to be 3-^/3 « 5.2 times larger than 
would be in principal necessary to extract the same infor- 
mation. Particularly, in the context of ab-initio molecu- 
lar dynamics (MD) [3 - 5 this results in a substantial addi- 
tional computational burden. The situation is even worse 
for a general orthorhomic unit cell, where the fraction of 



the exploited volume may become arbitrarily small. For 
quasi one-dimensional systems, such as nanowires or nan- 
otubes, this does indeed occur, just as in slab calculations 
of two-dimensional surfaces, or during shock compression 
simulations, where the unit cell is either rather prolate 
or oblate. 

For the simplest case of a cube in three dimensions, De- 
serno has shown how to extend the permissible radius to 
the maximum value given by the particular unit cell [6]. 
In this Brief Report, we will generalize this approach to 
arbitrary orthorhombic unit cells with lattice parameters 
(a, b,c). 

Let us start be rewriting the PCF in terms of 5- 
functions for the first n particles and integrating over 
all N particles: 
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The latter denotes the ensemble average of the quantity 
YYi = ]_5{Ri — R'j), using R[,...,R' N as integration vari- 
ables. Again, confining ourself to the most relevant case 



n = 2 and invoking isotropy, Eq. ( 2b I can be written as 
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where p = N/V is the particle density, while G(r) may be 
any arbitrary quantity for which to calculate the RDF. In 
a numerical calculation, the numerator G(r) is computed 
by averaging G(r) over shells of constant radius via a 
binning procedure, while the denominator A(r) plays the 
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role of a normalization and depends only on the geometry 
of the simulation box. However, due to the discretization 
of space, we cannot directly compute G(r), but solely its 
integrated value F(r) over the length of a bin Ar: 



dr'G(r') := F(r) 
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The choice of normalization for F(r) is arbitrary, except 
for the constraint that the normalization must reproduce 
Eq. ([3]), i.e. in the limit Ar — > it converges to the exact 
solution. The most straightforward choice of A(r) is 
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where V(r) :— dr' A(r') is the volume of the sphere 
that is intersecting with the simulation cell. Instead of 
Eq. A(r) x Ar can be employed as an alternative 
normalization. The fact that both are indeed genuine 
normalizations can be seen using the rule of l'Hopital: 
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I-VIII and Vj(r) through Vvin(r), respectively: 



/ 
II 
III 

IV 

V 
VI 
VII 
VIII 



< 2r < c 
c < 2r < b 

b < 2r < min(a, \J c 2 + b 2 ) 
min(a, \J c 2 + b 2 ) < 2r < max(a, \/c 2 + b 2 ) 
IV a :a< Vc 2 + b 2 , IV b : V ' + b 2 < a 
max(a, \/c 2 + b 2 ) <2r < \Jc 2 + a 2 



\J c 2 + a 2 < 2r < y/b 2 + a 2 



sjb 2 + a 2 < 2r < \/c 2 + b 2 + a 2 
\J c 2 +b 2 + a 2 < 2r 



The cases IV a and IVb are mutually exclusive, and it 
depends on the cell geometry, which one is appropriate. 
In cartesian coordinates, V(r) reads as 
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Since F(r) can be directly extracted from a MD or 
Monte Carlo simulation, in the following we will only 
discuss how to evaluate Eq. ([5]). To that extent, we will 
provide an explicit set of formulae to compute V(r) for 
the general case of an orthorhombic simulation cell. 

If the size of the simulation cell is infinite, the nor- 
malization function V(r) would simply correspond to 
the volume of a sphere, i.e. Vo(r) = ^7rr 3 . For a fi- 
nite simulation cell though, V(r) is the intersecting vol- 
ume of the sphere with the simulation box. This re- 
quires a distinction of cases for the different types of 
intersections, which complicates matters. The simplest 
solution is to circumvent this difficulty by confining r, 
so that < 2r < mm(L x , L y , L z ) holds. In this way 
Vo(r) = l^r 3 is still valid, but unfortunately, any addi- 
tional information outside the sphere is neglected. Thus, 
depending on the shape of the simulation cell, the uti- 
lized radius is at least -\/3 « 1.7 times smaller than ide- 
ally possible. However, the latter requires to compute 
V(r) for all different types of intersections. For a general 
orthorhombic simulation cell with dimensions a > b > c, 
we can distinguish between seven critical values of r: 

l<\<^<\Vc^b~ 2 <i^W 

< ^Vb 2 + a 2 < ^x/c 2 + b 2 +a 2 (7) 

This results in the following eight intervals for r and 
corresponding volume functions, which are denoted as 



where the integration limits are 

X(r) - min (r, |) , (9a) 
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respectively. 

Solving Eq. (JsJ) analytically for each of the aforemen- 
tioned cases, eventually leads to 
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FIG. 1. V(r) for an orthorhombic unit cell with dimensions 
2x3x5 (red) and 4x5x6 (blue). The dashed line indicates 
Vo(r) = §7rr 3 . 



FIG. 2. A(r) = d r V(r) for an orthorhombic unit cell with 
dimensions 2x3x5 (red) and 4x5x6 (blue). The dashed 
line indicates Ao(r) — Airr 2 . 



where we have introduced the shorthand notations 



V K {a) := it 
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and 
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Vn(P, 7) := g [H(P, 7) + #(7, £)] , with (10b) 
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Herein, Vk{o) denotes the volume of a spheric dome of 
radius one that is cutoff at a height z — a, whereas 
Vjv(/3, 7) corresponds to the volume of a wedge cut from 
the unit sphere, whose planar sides lie at x = and 
y = 7. Since the auxiliary functions Vk(oc) and H{uj,t) 
are only needed for a € [0, 1] and oj 2 +t 2 € [0, 1], respec- 
tively, it is possible to exploit this in the interest of an 
efficient implementation. 

In Fig. [l] and [2] the eventual normalization function 
V(r) and the denominator A(r) = d r V{r) of Eq. ^ 
are shown for two sample simulation cells and com- 
pared to Vo(r) — |7rr 3 and A (r) = Airr 2 , respec- 
tively. As can be directly seen, the accuracy of the 
functions Vo{r) and A$(r) immediately deteriorates for 
r > min(L a; , L y , L z )/2. Although, using the present com- 
putational method, it is possible to exploit the whole 
volume up to the maximum permissible radius r max = 
V a 2 + b 2 + c 2 /2, the statistics will become inferior with 



increasing radius. Moreover, for r — > r max the evaluation 
of Eq. ([3| will not only be affected by statistical uncer- 
tainties, but also by a small bias due to the fact that in 
this case the normalization function is vanishing. 

For the purpose of demonstrating our new calcula- 
tion method, we present here shock compression simu- 
lations of liquid hydrogen to illustrate that our method 
works well for highly prolate and oblate cell geome- 
tries. However, instead of simulating a planar shock 
wave within a large computational box with many atoms, 
we have employed the multiscale shock-wave simulation 
technique (MSST) [5], which is based on MD as well 
as the Navier-Stokes equations for compressible flow. In 
this way, the simulation cell follows a Lagrangian point 
through the shock wave as if the shock were passing over 
it, instead of directly simulating the shock- wave itself. 
This allows to perform shock-compression simulations 
with a sufficiently small number of atoms, so that it be- 
came feasible to perform MSST simulations in conjunc- 
tion with ab-initio MD. 

Due to the fact that a direct Born-Oppenheimer MD 
(BOMD) simulation, where the total energy functional 
is fully optimized in every MD step, would have been 
prohibitive, we employ here the recently devised "Car- 
Parrinello-like approach to Born-Oppenheimer MD" of 
Kiihne et al. [1], which has already demonstrated its su- 
perior efficiency [HI HO] • In the spirit of the original Car- 
Parrinello MD (CPMD) method |3] during the dynamics, 
the electronic wave functions are not self-consistently op- 
timized. Nevertheless, at variance to CPMD, in this ap- 
proach the fictitious Newtonian dynamics of the electrons 
and ions is substituted by a similarly coupled electron-ion 
MD, which does not require the definition of a fictitious 
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mass parameter, but at the same time keeps the electrons 
very close to their instantaneous electronic ground state. 
As a consequence, the time step can be chosen up to 
the ionic resonance limit, while simultaneously preserv- 
ing the efficiency of CPMD. As a consequence, the best 
aspects of the BOMD and CPMD schemes are unified, 
which not only extends the scope of either method but 
allows for ab-initio simulations previously thought to be 
not feasible. 

The MSST simulation has been performed within den- 
sity functional theory using the mixed Gaussian and 
plane wave [11] code CP2K/Quickstep [12], where the 
electron density is represented by a plane wave basis set, 
while the orbitals are expanded in Gaussians. Together 
with efficient transformation methods to switch between 
one representation or the other and advanced multigrid, 
screening as well as sparse matrix methods, an efficient 
linear-scaling evaluation of the Hamilton matrix is ob- 
tained. Efforts towards a full linear-scaling algorithm 
are underway [13] • Here the orbitals are described by an 
accurate double-^ set with one set of polarization func- 
tions (DZVP) [13], whereas a density cutoff of 200 Ry is 
employed for the charge density. The unknown exchange 
and correlation potential is substituted by the PBE gen- 
eralized gradient approximation 15j. The interactions 
between the valence electrons and the ionic cores are de- 
scribed by rather hard norm-conserving pseudopotentials 

dEinzi. 

The system consists of 768 hydrogen atoms in an or- 
thorombic box of initial dimensions 68.2 A x 11.0 A x 
11.0 A. Before starting the actual MSST simulation, the 
sample has been equilibrated in the NPT ensemble at a 
temperature of 100 K and a pressure of 1 GPa. Since we 
are dealing with a rather large disordered system at finite 
temperature, the Brillouin zone is sampled at the T-point 
only. The MSST equations of motion for the nuclei and 
volume of the simulation cell are integrated using a dis- 
cretized time step of 0.1 fs and a cell mass of 7 x 10 7 a.u. 
(= mass x length -4 ) to constrain the stress in the shock 
propagation direction to the Rayleigh line and the energy 
of the system to the Hugonoit energy condition. 

Upon applying a shock velocity of 35 km/s along the 
x-axis, the system undergoes a 9.2-fold compression in 
volume along the same direction resulting in a simula- 
tion box of dimensions 7.4 Ax 11.0 Axll.O A. In doing 
so the pressure increases till 650 GPa, while the tempera- 
ture raises up to 8500 K. The according RDF's before and 
after the shock are shown in Fig. [3] and |1J which can be 
compared to the conventional <?(r) in the inset. Using the 
present scheme, it is possible to extract as much informa- 
tion that otherwise would require a simulation cell which 
is larger by a factor of 38 and 1.5, respectively. In Fig. [3] 
it is also also observed that liquid hydrogen undergoes 
a molecular-atomic transition to a metal, which is con- 
sistent with the findings of Nellis and coworkers [T5J [TH] • 
However, further details and additional simulations will 
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FIG. 3. The RDF of liquid molecular Hydrogen before the 
shock-compression at a temperature of 100 K and a pressure 
of 1 GPa. 
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FIG. 4. The RDF of liquid atomic Hydrogen after the shock- 
compression at a temperature of 8500 K and a pressure of 
650 GPa. 



be published elsewhere. 

We conclude by noting that the observed increase in 
terms of a accessible volume does neither mitigate single 
particle finite size effects, nor does it substitute a proper 
finite size scaling [20] . 
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